

// /// explicit specialization of the whole class for complex
// // template <> struct sgArpackPrivate<Complex_t>;
// template <>
// class ALGEBRA_EXPORT sgArpack<Complex_t>
// {
//   public:
//     sgArpack ();
//     ~sgArpack ();

//     struct Option
//     {
//         /*
//         c  WHICH   Character*2.  (INPUT)
//         c          Specify which of the Ritz values of OP to compute.
//         c
//         c          'LA' - compute the NEV largest (algebraic) eigenvalues.
//         c          'SA' - compute the NEV smallest (algebraic) eigenvalues.
//         c          'LM' - compute the NEV largest (in magnitude) eigenvalues.
//         c          'SM' - compute the NEV smallest (in magnitude) eigenvalues.
//         c          'BE' - compute NEV eigenvalues, half from each end of the
//         c                 spectrum.  When NEV is odd, compute one more from the
//         c                 high end than from the low end.
//         */
//         std::string which = "SA";  // Specify which of the Ritz values of OP to compute.
//         int         nev   = 1;     // Number of eigenvalues of OP to be computed. 0 < NEV < N.
//         Complex_t   sigma = {0.0, 0.0};   // shifts
//     };

//     /// the default option to control the eigenvalue solver.
//     Option getDefaultOption ();

//     /// Number of eigenvalues of OP to be computed. 0 < NEV < N.
//     void setOption (const Option& opt);

//     /**
//      * @brief solve the generalized eigenvalue problem A*x = lambda*B*x.
//      *
//      * Mode 3:  K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite.
//      *          ===> OP = (inv[K - sigma*M])*M  and  B = M.
//      *          ===> Shift-and-Invert mode
//      * @return [eigenvectors, eigenvalues]:
//      *            eigenvectors = V(:. 0:nev),
//      *            eigenvalues = V(0:nev, nev)
//      */
//     ComplexMat_t solve (const ComplexSpMat_t& A, const ComplexSpMat_t& B);

//     // /// A matlab-like signature to solve the generalized eigenvalue problem
//     // // see http://www.ece.northwestern.edu/local-apps/matlabhelp/techdoc/ref/eig.html
//     // void solve (const SpMat_t& A, const SpMat_t& B, Mat_t& V, Mat_t& D);

//     /// report whether previous computation was successful.
//     Info_t info () const;

//   private:
//     struct sgArpackPrivate<Complex_t>* imp = nullptr;
// };

// // Give the explicit instantiatioon to hide the underlying implementation.
// // The compiler will use the specialization version of Complex_t
// template class ALGEBRA_EXPORT sgArpack<Complex_t>;

// sgArpack<Complex_t>::sgArpack ()
// {}

// sgArpack<Complex_t>::~sgArpack ()
// {}

// sgArpack<Complex_t>::Option sgArpack<Complex_t>::getDefaultOption ()
// {
//     static sgArpack::Option opt;

//     return opt;
// }

// void sgArpack<Complex_t>::setOption (const sgArpack<Complex_t>::Option& opt)
// {}

// ComplexMat_t sgArpack<Complex_t>::solve (const ComplexSpMat_t& A, const ComplexSpMat_t& B)
// {
//     ComplexMat_t E;

//     return E;
// }

// Info_t sgArpack<Complex_t>::info () const
// {
//     return ( Info_t ) imp->info;
// }

// TEST_F (ODETest, NewmarkTest)
// {
//     typedef Eigen::Triplet<Real_t> Triplet;

//     // ODE System
//     ForcedDynamicSystem  ode;
//     Real_t               m       = 1.0;
//     Real_t c       = 1.0;
//     Real_t k       = 1.0;
//     Real_t F0      = 1.0;
//     Real_t omega   = 1.0;

//     int     n = 1;  // the number of equations in the system
//     SpMat_t M (n, n), C (n, n), K (n, n);
//     Vec_t   F (n), U (n), V (n), A (n);

//     //-- [M]
//     {
//         std::vector<Triplet> triplets;

//         triplets.push_back (Triplet (0, 0, m));
//         M.setFromTriplets (triplets.begin (), triplets.end ());
//     }

//     //-- [C]
//     {
//         std::vector<Triplet> triplets;
//         triplets.push_back (Triplet (0, 0, c));
//         C.setFromTriplets (triplets.begin (), triplets.end ());
//     }

//     //-- [K]
//     {
//         std::vector<Triplet> triplets;
//         triplets.push_back (Triplet (0, 0, m));
//         K.setFromTriplets (triplets.begin (), triplets.end ());
//     }

//     // initial conditions
//     Real_t u0 = 1.0;
//     Real_t v0 = 0.0;
//     U[0] = u0;
//     V[0] = v0;

//     // ODE Solver

//     Newmark solver;
//     solver.Initialize (ode);
//     int     steps = 100;
//     Real_t  t0    = 0.0;  // start time
//     Real_t  tf    = 1.0;  // final time
//     Real_t  dt    = (tf -t0)/ steps;  // time step

//     // numerical solution
//     std::vector<double> numerical_solution;

//     numerical_solution.push_back (u0);
//     for (int i = 1; i < steps; ++i)
//     {
//         Real_t t = i * dt;

//         solver.SetTime (t);
//         solver.SetTimeStep (dt);

//         F[0] = F0 * sin (omega * t);

//         solver.Step (M, C, K, F, U, V, A);

//         numerical_solution.push_back (U[0]);
//     }

//     // analytical solution
//     Real_t              omega_n = sqrt (k / m);
//     Real_t              zeta    = c / (2.0 * m * omega_n);
//     Real_t              zeta2   = zeta * zeta;
//     Real_t              omega_d = omega_n * sqrt (1.0 - zeta2);
//     Real_t              gama    = omega / omega_n;
//     Real_t              gama2   = gama * gama;
//     Real_t              X       = F0 / k * sqrt ((1.0 - gama2) * (1.0 - gama2) + 4 * zeta2 * gama2);
//     Real_t              alpha   = atan (2.0 * zeta * gama / (1.0 - gama2));

//     std::vector<double> analytical_solution;
//     for (int i = 0; i < steps; ++i)
//     {
//         Real_t t     = i * dt;
//         Real_t theta   = omega * t;
//         Real_t theta_n = omega_n * t;
//         Real_t theta_d = omega_d * t;
//         Real_t u;

//         // term1
//         Real_t fact1 = (v0 + zeta * omega_n * u0);
//         fact1        = fact1 * sin (theta_d) / omega_d;
//         fact1        = fact1 + u0 * cos (theta_d);
//         u            = exp (-zeta * theta_n) * fact1;

//         // term2
//         Real_t fact2 = zeta * sin (alpha) - omega * cos (alpha);
//         fact2        = fact2 * sin (theta_d) / omega_d;
//         fact2        = fact2 + sin (alpha) * cos (theta_d);
//         u            = u + X * exp (-zeta * theta_n) * fact2;

//         // term3
//         u = u + X * sin (theta - alpha);

//         analytical_solution.push_back (u);
//     }

// #ifdef SG_VERBOSE
//     std::cout << "numerical/analytical solution = " << std::endl;
//     for (int i = 0; i < steps; ++i)
//     {
//         std::cout << numerical_solution[i] << ", " << analytical_solution[i] << std::endl;
//     }
// #endif

//     // testing the tolerance between the numerical/analytical solution
//     for (int i = 0; i < steps; ++ i)
//     {
//         ASSERT_NEAR (numerical_solution[i], analytical_solution[i], TOLERANCE);
//     }
// }

#ifdef SG_USE_MKL
#include "mkl.h"
//
// Some references about oneMKL,
//
// 1. The matrix storage in Eigen
// The SparseMatrix and SparseVector classes take three template arguments:
// the scalar type (e.g., double) the storage order (ColMajor or RowMajor,
// the default is ColMajor) the inner index type (default is int).
// see https://eigen.tuxfamily.org/dox/group__TutorialSparse.html
//
// 2. Computes matrix-vector product for a sparse matrix in the CSC format (deprecated).
// void mkl_dcscmv (const char    *transa,
//                  const MKL_INT *m,
//                  const MKL_INT *k,
//                  const double  *alpha,
//                  const char    *matdescra,
//                  const double  *val,
//                  const MKL_INT *indx,
//                  const MKL_INT *pntrb,
//                  const MKL_INT *pntre,
//                  const double  *x,
//                  const double  *beta,
//                  double        *y)
// see https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-2/mkl-cscmv.html
//
// 3. Creates a handle for a CSC format matrix.
// sparse_status_t mkl_sparse_d_create_csc (sparse_matrix_t * A,
//                                          const sparse_index_base_t indexing,
//                                          const MKL_INT             rows,
//                                          const MKL_INT             cols,
//                                          MKL_INT                  *cols_start,
//                                          MKL_INT                  *cols_end,
//                                          MKL_INT                  *row_indx,
//                                          double                   *values)
// see https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-2/mkl-sparse-create-csc.html
//
// 4. Computes a sparse matrix- vector product.
// sparse_status_t mkl_sparse_d_mv (const sparse_operation_t  operation,
//                                  const double              alpha,
//                                  const sparse_matrix_t     A,
//                                  const struct matrix_descr descr,
//                                  const double*             x,
//                                  const double              beta,
//                                  double*                   y)
// see https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-2/mkl-sparse-mv.html
//

// Info_t SparseBlas::dot (const Vec_t& X, const Vec_t& Y, Real_t& R)
// {
//     int exit_status = 0;

//     MKL_INT n = X.size ();

//     const double* x = X.data ();
//     const double* y = Y.data ();

//     MKL_INT incx = 1, incy = 1;

//     R = cblas_ddot (n, x, incx, y, incy);

//     return ( Info_t ) exit_status;
// }

// Info_t dot (const ComplexVec_t& X, const ComplexVec_t& Y, Complex_t& R)
// {
//     int exit_status = 0;

//     MKL_INT n = X.size ();

//     auto  x    = X.data ();
//     auto  y    = Y.data ();
//     void* dotu = &R;

//     MKL_INT incx = 1, incy = 1;

//     cblas_zdotu_sub (n, x, incx, y, incy, dotu);

//     return ( Info_t ) exit_status;
// }

Info_t SparseBlas::multiply (const SpMat_t& A, Vec_t& X, Vec_t& Y)
{
#if 1
    int exit_status = 0;
    // Specifies the operation.
    // If      transa = 'N' or 'n',
    //   then y :    = alpha *A *x + beta *y
    // If transa = 'T' or 't' or 'C' or 'c',
    //   then y : = alpha * AT * x + beta * y
    const char* transa       = "N";
    MKL_INT     nnz          = A.nonZeros ();  // the number of non zero coefficients
    MKL_INT     m            = A.rows ();      // Number of rows of the matrix A
    MKL_INT     k            = A.cols ();      // Number of columns of the matrix A
    double      alpha        = 1.0;            // Specifies the scalar alpha.
    char        matdescra[6] = { '\0' };       // Array of six elements, specifies properties of the matrix used for operation.
    // Array of six elements, specifies properties of the matrix used for operation. Only first four array elements are used,
    // their possible values are given in Table “Possible Values of the Parameter matdescra (descra)”. see
    // https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-2/interface-consideration.html#TBL2-6
    matdescra[0] = 'G';                 // matrix structure: general
    matdescra[1] = 'L';                 // upper/lower triangular indicator: ignored
    matdescra[2] = 'N';                 // main diagonal type: ignored
    matdescra[3] = 'C';                 // type of indexing: zero-based indexing

    const double* val = A.valuePtr ();  // Array containing non-zero elements of the matrix A.

    // For ILP64 interface, MKL_INT will be defined as __int64 or long long int.
    // Which is different with SpMat_t::Index. Thus, resort to deep-copy for compatibilities.

    //
    // The indx/pntrb/pntre description on the official website seems wrong.
    // For zero-based indexing, array containing the row indices for each non-zero element of the matrix A.
    // see https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-2/mkl-cscmv.html
    // see https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-2/mkl-csrmv.html
    //
    //

    // indx
    // For one-based indexing, array containing the column indices plus one for each non-zero element of the matrix A.
    // For zero-based indexing, array containing the column indices for each non-zero element of the matrix A.
    // Its length is equal to length of the val array.
    std::vector<MKL_INT> indx (nnz);
    for (auto i = 0; i < nnz; ++i)
    {
        indx[i] = A.innerIndexPtr ()[i];
    }

    // pntrb
    // Array of length k.
    // This array contains column indices, such that pntrb[i] - pntrb[0] + 1 is the first index of column i in the arrays val and indx.
    std::vector<MKL_INT> pntrb (k);
    for (auto i = 0; i < k; ++i)
    {
        pntrb[i] = A.outerIndexPtr ()[i];
    }

    // pntre
    // Array of length k.
    // For one-based indexing this array contains column indices, such that pntre[i] - pntrb[1] is the last index of column i in the arrays
    // val and indx. For zero-based indexing this array contains column indices, such that pntre[i] - pntrb[1] - 1 is the last index of
    // column i in the arrays val and indx.
    std::vector<MKL_INT> pntre (k);
    for (auto i = 0; i < k; ++i)
    {
        pntre[i] = A.outerIndexPtr ()[i + 1];
    }

    const double* x = X.data ();  // Array, size at least k if transa = 'N' or 'n' and at least m otherwise. On entry, the array x must
                                  // contain the vector x.

    double beta = 0.0;            // Specifies the scalar beta.

    Y.resize (m);
    double* y = Y.data ();        // Array, size at least m if transa = 'N' or 'n' and at least k otherwise. On entry, the array y must
                                  // contain the vector y.

    // call the oneMKL Sprase BLAS routine
    mkl_dcscmv (transa, &m, &k, &alpha, matdescra, val, &indx[0], &pntrb[0], &pntre[0], x, &beta, y);

    return ( Info_t ) exit_status;
#else
    // snippets from Intel /opt/oneAPI/mkl/2024.2/share/doc/mkl/examples/c/sparse_blas/source/sparse_convert.c
    int exit_status = 0;

    sparse_status_t status;

    // Structure with sparse matrix stored in CSR format
    sparse_matrix_t cscA;

    sparse_index_base_t indexing = SPARSE_INDEX_BASE_ZERO;  // Indicates how input arrays are indexed:
                                                            // Zero-based (C-style) indexing: indices start at 0.
    MKL_INT nnz = A.nonZeros ();                            // the number of non zero coefficients
    MKL_INT m   = A.rows ();                                // Number of rows of the matrix A
    MKL_INT k   = A.cols ();                                // Number of columns of the matrix A

    // For ILP64 interface, MKL_INT will be defined as __int64 or long long int.
    // Which is different with SpMat_t::Index. Thus, resort to deep-copy for compatibilities.

    // cols_start
    // Array of length at least m. This array contains col indices, such that cols_start[i] - ind is the first index of
    // col i in the arrays values and row_indx. ind takes 0 for zero-based indexing and 1 for one-based indexing.
    std::vector<MKL_INT> cols_start (k);
    for (auto i = 0; i < k; ++i)
    {
        cols_start[i] = A.outerIndexPtr ()[i];
    }

    // cols_end
    // Array of at least length m. This array contains col indices, such that cols_end[i] - ind - 1 is the last index of
    // col i in the arrays values and row_indx. ind takes 0 for zero-based indexing and 1 for one-based indexing.
    std::vector<MKL_INT> cols_end (k);
    for (auto i = 0; i < k; ++i)
    {
        cols_end[i] = A.outerIndexPtr ()[i + 1];
    }

    // row_indx
    // For one-based indexing, array containing the row indices plus one for each non-zero element of the matrix A.
    // For zero-based indexing, array containing the row indices for each non-zero element of the matrix A.
    // Its length is at least cols_end[cols - 1] - ind. ind takes 0 for zero-based indexing and 1 for one-based indexing.
    std::vector<MKL_INT> row_indx (nnz);
    for (auto i = 0; i < nnz; ++i)
    {
        row_indx[i] = A.innerIndexPtr ()[i];
    }

    // values
    // Array containing non- zero elements of the matrix A.Its length is equal to length of the row_indx array.
    double* values = const_cast<SpMat_t&> (A).valuePtr ();  // Array containing non-zero elements of the matrix A.

    // x
    // Array of size equal to the number of columns, k of A if operation = SPARSE_OPERATION_NON_TRANSPOSE and at least
    // the number of rows, m, of A otherwise.On entry, the array must contain the vector x.
    const double* x = X.data ();

    // beta
    // Specifies the scalar beta.
    double beta = 0.0;

    // y
    // Array with size at least m if operation = SPARSE_OPERATION_NON_TRANSPOSE and at least k otherwise.On entry,
    // the array y must contain the vector y.Array of size equal to the number of rows,
    // m of A if operation = SPARSE_OPERATION_NON_TRANSPOSE and at least the number of columns, k,
    // of A otherwise.On entry, the array y must contain the vector y.
    Y.resize (m);
    double* y = Y.data ();

    sparse_operation_t  operation = SPARSE_OPERATION_NON_TRANSPOSE;  // Specifies operation op() on input matrix.
    double              alpha     = 1.0;                             // Specifies the scalar alpha.
    struct matrix_descr descrA;                                      // Structure specifying sparse matrix properties
    descrA.type = SPARSE_MATRIX_TYPE_GENERAL;                        // The matrix is processed as is.

    // The mkl_sparse_ ? _create_csc routine creates a handle for an m - by - k matrix A in CSC format.
    status = mkl_sparse_d_create_csc (&cscA, indexing, m, k, &cols_start[0], &cols_end[0], &row_indx[0], values);

    if (status != SPARSE_STATUS_SUCCESS)
    {
        printf (" Error in mkl_sparse_d_create_csc: %d \n", status);
        exit_status = 1;
        goto EXIT;
    }

    status = mkl_sparse_d_mv (operation, alpha, cscA, descrA, x, beta, y);
    if (status != SPARSE_STATUS_SUCCESS)
    {
        printf (" Error in mkl_sparse_d_mv: %d \n", status);
        exit_status = 1;
        goto EXIT;
    }

EXIT:
    return ( Info_t ) exit_status;
#endif
}

Info_t SparseBlas::multiply (const ComplexSpMat_t& A, ComplexVec_t& X, ComplexVec_t& Y)
{
    int exit_status = 0;
    // Specifies the operation.
    // If      transa = 'N' or 'n',
    //   then y :    = alpha *A *x + beta *y
    // If transa = 'T' or 't' or 'C' or 'c',
    //   then y : = alpha * AT * x + beta * y
    const char*   transa       = "N";
    MKL_INT       nnz          = A.nonZeros ();  // the number of non zero coefficients
    MKL_INT       m            = A.rows ();      // Number of rows of the matrix A
    MKL_INT       k            = A.cols ();      // Number of columns of the matrix A
    MKL_Complex16 alpha        = { 1.0, 0.0 };   // Specifies the scalar alpha.
    char          matdescra[6] = { '\0' };       // Array of six elements, specifies properties of the matrix used for operation.
    // Array of six elements, specifies properties of the matrix used for operation. Only first four array elements are used,
    // their possible values are given in Table “Possible Values of the Parameter matdescra (descra)”. see
    // https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-2/interface-consideration.html#TBL2-6
    matdescra[0] = 'G';                                           // matrix structure: general
    matdescra[1] = 'L';                                           // upper/lower triangular indicator: ignored
    matdescra[2] = 'N';                                           // main diagonal type: ignored
    matdescra[3] = 'C';                                           // type of indexing: zero-based indexing

    const MKL_Complex16* val = ( MKL_Complex16* ) A.valuePtr ();  // Array containing non-zero elements of the matrix A.

    // For ILP64 interface, MKL_INT will be defined as __int64 or long long int.
    // Which is different with SpMat_t::Index. Thus, resort to deep-copy for compatibilities.

    //
    // The indx/pntrb/pntre description on the official website seems wrong.
    // For zero-based indexing, array containing the row indices for each non-zero element of the matrix A.
    // see https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-2/mkl-cscmv.html
    // see https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-2/mkl-csrmv.html
    //
    //

    // indx
    // For one-based indexing, array containing the column indices plus one for each non-zero element of the matrix A.
    // For zero-based indexing, array containing the column indices for each non-zero element of the matrix A.
    // Its length is equal to length of the val array.
    std::vector<MKL_INT> indx (nnz);
    for (auto i = 0; i < nnz; ++i)
    {
        indx[i] = A.innerIndexPtr ()[i];
    }

    // pntrb
    // Array of length k.
    // This array contains column indices, such that pntrb[i] - pntrb[0] + 1 is the first index of column i in the arrays val and indx.
    std::vector<MKL_INT> pntrb (k);
    for (auto i = 0; i < k; ++i)
    {
        pntrb[i] = A.outerIndexPtr ()[i];
    }

    // pntre
    // Array of length k.
    // For one-based indexing this array contains column indices, such that pntre[i] - pntrb[1] is the last index of column i in the arrays
    // val and indx. For zero-based indexing this array contains column indices, such that pntre[i] - pntrb[1] - 1 is the last index of
    // column i in the arrays val and indx.
    std::vector<MKL_INT> pntre (k);
    for (auto i = 0; i < k; ++i)
    {
        pntre[i] = A.outerIndexPtr ()[i + 1];
    }

    const MKL_Complex16* x = ( MKL_Complex16* ) X.data ();  // Array, size at least k if transa = 'N' or 'n' and at least m otherwise. On
                                                            // entry, the array x must contain the vector x.

    MKL_Complex16 beta = { 0.0, 0.0 };                      // Specifies the scalar beta.

    Y.resize (m);
    MKL_Complex16* y = ( MKL_Complex16* ) Y.data ();  // Array, size at least m if transa = 'N' or 'n' and at least k otherwise. On entry,
                                                      // the array y must contain the vector y.

    // call the oneMKL Sprase BLAS routine
    mkl_zcscmv (transa, &m, &k, &alpha, matdescra, val, &indx[0], &pntrb[0], &pntre[0], x, &beta, y);

    return ( Info_t ) exit_status;
}

#else
// Info_t SparseBlas::dot (const Vec_t& X, const Vec_t& Y, Real_t& R)
// {
//     int exit_status = 0;

//     R = X.dot (Y);

//     return ( Info_t ) exit_status;
// }

// Info_t dot (const ComplexVec_t& X, const ComplexVec_t& Y, Complex_t& R)
// {
//     int exit_status = 0;

//     R = X.dot (Y);

//     return ( Info_t ) exit_status;
// }

Info_t SparseBlas::multiply (const SpMat_t& A, Vec_t& X, Vec_t& Y)
{
    Y = A * X;

    return Info_t (0);
}

Info_t SparseBlas::multiply (const ComplexSpMat_t& A, ComplexVec_t& X, ComplexVec_t& Y)
{
    Y = A * X;

    return Info_t (0);
}

#endif